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We present two experimental schemes that can be used to implement the Factorized Quantum 
Lattice-Gas Algorithm for the ID Diffusion Equation with Persistent-Current Qubits. One scheme 
involves biasing the PC Qubit at multiple flux bias points throughout the course of the algorithm. 
An implementation analogous to that done in Nuclear Magnetic Resonance Quantum Computing 
is also developed. Errors due to a few key approximations utilized are discussed and differences 
between the PC Qubit and NMR systems are highlighted. 
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INTRODUCTION 



o 

. Most algorithms designed for quantum computers will not best their classical counterparts until they are imple- 
mented with thousands of qubits. For example, the factoring of binary numbers with a quantum computer is estimated 
to be faster than a classical computer only when the length of the number is greater than about 500 digits 1 . Account- 
ing for error correction circuitry^ would bring the size of the needed quantum computer to be in the thousands of 
qubits. In contrast, the Factorized Quantum Lattice-Gas Algorithm (FQLGA)^ for fluid dynamics simulation, even 
when run on a quantum computer significantly smaller than the one just discussed, has significant advantages over 
& \ its classical counterparts. 

The FQLGA is the quantum version of classical lattice-gases (CLG)^. CLG are an extension of classical cellular 
automata with the goal of simulating fluid dynamics without reference to specific microscopic interactions. The 
binary nature of the CLG lattice variables is replaced for the FQLGA by the Hilbert space of a two-level quantum 
system. The results of this replacement are similar to that of the lattice-Boltzmann model, but with a few significant 
differences^. The first is the exponential decrease in required memory. The second is the ability to simulate arbitrarily 
■ small viscosities. 

As of today there is a plethora of qubits to choose from when designing a quantum computer, and a promising 
class is superconducting qubits based on Josepshon junction circuits&L&Siifi. One major advantage of any of these 
superconducting systems is the ability to precisely engineer the quantum Hamiltonian, which extends from single 
qubit design to multi-qubit coupling arrangements to measurement engineering. The quantum computer considered 
here will be built using the Persistent-Current Qubit (PC Qubit )£. 

The goal of this paper is to show how one can implement a one dimensional version of the FQLGA with the PC Qubit. 
To this end we will begin by reviewing the algorithm, specifically the one that simulates the diffusion equation, without 
a loss of generality in understanding the essence of the algorithm or its general requirements. We will then review the 
PC qubit and show explicitly how to implement the algorithm with this system. Some important differences between 
the PC qubit and the two-state system studied in Nuclear Magnetic Resonance Quantum Computation (NMRQC)ii 
will be shown to allow for some interesting new techniques in implementing quantum logic. We will also show how 
to implement the algorithm with the PC qubit in a very analogous way to NMRQC schemes^, with a few significant 
differences. 



II. FQLGA FOR THE ID DIFFUSION EQUATION 



The first thing one must do in the FQLGA is to define a lattice. Each lattice point n will represent a unique position 
in the simulated fluid. The simulation will contain a finite number of lattice points, hence space is discretized in the 
simulation. 
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Next one must encode the mass density p of the fluid at each lattice site. In the FQLGA this is done by building at 
each lattice site a set {i} of coupled qubits. Each qubit represents the motion of particles on the microscopic level in 
one of a finite set of directions. For the diffusion equation in one dimension, at any point in your fluid, there are only 
two possible directions for each particle to be moving, to the left and to the right. Hence, only two qubits are needed 
to specify the mass density p n at each lattice site. This intuitive reasoning does not extrapolate to higher dimensional 
simulations because even in two dimensions there would be an infinite number of directions particles could travel in. 
In higher dimensions one must adhere to much more mathematical conditions to decide on the small set of directions 
one must include for a faithful simulation^. The probability P of a particle to be participating in the motion assigned 
to each qubit will be encoded in the probability amplitude of the qubit being in its excited state The state of a 
qubit is thus set to 



where i is the qubit index, n is the lattice site index, and |0) is the ground state of the qubit. For the one dimensional 
problem considered here, i — {1,2} and n — {1, N} where N is the number of lattice sites used in the simulation. One 
can easily conceive of fluids of multiple phases with multiple types of interactions even in one dimension, in which 
the size of {i} would be much larger, but this will not be considered here. The mass density p is then calculated 
by summing the occupation probabilities for all qubits at a node. At time t=0 in a ID simulation the occupation 
probabilities P\ n and P^ 1 are set to p n /2, which is the condition for local equilibrium in the fluidi^. 

Now that the fluid is initialized, one must account for the interaction of particles in the fluid. These collisions 
are encoded by the application of a unitary transformation to the coupled systems at each lattice site. For the ID 
diffusion equation this unitary transformation is 



The basis for computation is the set of four product states: |0)|0), representing no particles at the site, |0)|1), 
representing the existence of only a particle moving to the right at the site, |1)|0), representing the existence of only 
a particle moving to the left at the site, and |1)|1), representing particles moving in both directions at the site. To 
conserve particle number there can be coupling only between the middle two states. The identity transformation on 
the first and last states corresponds to no collisions and a perfectly elastic collision respectively. Transformation of 
the middle two states was something that never existed in the classical algorithm because there was no superposition 
of these two states. 

After collision the states of the qubits at each lattice site are in general entangled, and we denote that state as 
|T"). The state of each qubit is then measured, denoted by |xi"), and the process described thus far is repeated many 
times to achieve an ensemble average. Upon completion of these measurements one will have found the post-collision 
outgoing occupation probabilities, denoted by once again. Note that the occupation probabilities now represent 
something very different than before the collision. The particles have now interacted and are ready to move to the 
next lattice site. 

One must now "stream" the occupation probabilities to their new lattice sites. This is done in a classical computer 
by storing the occupation probabilities at each lattice site that are coming from adjacent lattice sites due to colli- 
sions. More precisely, P\ n becomes P\ n+l and P2 1 becomes iV* - • Periodic boundary conditions are assumed when 
streaming at the edges of the fluid. 

To find the mass density p n at t=l one simply adds the occupation probability for both qubits at site n once 
streaming has been done. One time step of the algorithm has now been completed. To simulate the next time 
step simply start the above procedure all over again except now setting the initial states with the new occupation 
probabilities just found. 

The algorithm can be summarized by four major steps, which are illustrated in figure ^ The first step encodes 
the initial state of the fluid by quantum mechanically setting the state l^™) of each qubit at each lattice site. The 
second step transforms the two-qubit product state at each lattice site to in general an entangled state, whose state is 
denoted by |T"). Third one makes a projective measurement of the post-collision states |xi"), and one must repeat 
the first three steps to find the outgoing occupation probabilities P/\ In the fourth and final step one streams the 
mass density with the appropriate post-collision occupation probabilities, from the left with particles representing 
positive momentum, and from the right with particles representing negative momentum, and the mass density is 
calculated. Subsequent time steps are identical except for a change in the initial mass density profile, i.e., initial qubit 
states in the first step. 




(1) 




(2) 
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FIG. 1: General summary of the four major steps that comprise one time step of the ID FQLGA fluid dynamics simulation. 
The sequence of initialization of mass density, collision of particles, and measurement of post-collision states is repeated many 
times to make an ensemble measurement. Propagation between collisions is accomplished by storing the adjacent occupation 
probabilities for a given site in a classical computer, where the mass density is then calculated for this time step. Subsequent 
time steps utilize these "streamed" occupations when initializing again for the next set of collisions and "streaming" . 



III. PERSISTENT-CURRENT QUBIT 

The fundamental unit of quantum logic we will use to implement the algorithm is the PC Qubhiii It consists of a 
superconducting loop that is interrupted by three Josephson junctions, pictured as x's in figure Ufa). The magnetic 
flux <f> is the only control field for our qubit, and as shown in the figure, is usually denoted by / = <f>/$ where 
$ = h/2e is a single flux quantum, h is Planck's constant, and e is the magnitude of the charge of an electron. 
Physically, a Josephson junction is a small layer of insulator sandwiched between superconductors, so our system is 
a superconducting loop interrupted by three layers of insulator about lnm thick. For single qubit manipulation the 
magnetic flux through the loop will be modified. The flux seen by a DC SQUID magnetometer, a combination of 
applied flux and qubit-induced flux, will serve as our measurement variable. 

The Hamiltonian of the qubit is derived by considering a circuit element model of our system, which consists of 
three Josepshon junctions, where two junctions have the same cross-sectional area, and the third is smaller by a factor 
of a. The constituent relations for an ideal Josephson junction are 

I = I c sin(<^) (3a) 

where I is the current through the junction, V is the voltage across the junction, I c is the maximum current the 
junction can hold without a voltage appearing across it, (p=0i-02, and 6% t 2 is the phase of the plane wave macroscopic 
wavefunction that characterizes the superconductor condensate on the +,- side of the junction respectively. Note that 
I c is a function linear in the cross-sectional area of the junction, and hence the third junction has a lower I c by a 
factor of a. 

The energy associated with an ideal Josephson junction is found by integrating the power from time t = to some 
final time t a , which is equivalent to an integral from zero phase to some phase (p. The energy it takes to set the phase 
of a Josephson junction to <p is 

E = f ° (I e wa <f /)(p.^L)dt = %^ [* smip'dtf = Ej(l - cose/?) (4) 



' 2?r dt ' 2tt 



where Ej = $ I c /2ir. 
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FIG. 2: (a) Schematic drawing of the PC Qubit. The x's represent Josephson junctions, with all connecting leads made of 
the same superconductor that is part of the junctions. The sign conventions chosen when summing phases are shown, and the 
magnetic flux penetrating the loop (in units of $0) is labeled by /. (b) The potential energy of the full Hamiltonian for the 
PC Qubit is plotted when the system is biased at / = 0.495$ o - The phase particle sees an infinite 2D lattice with unit cells 
resembling a double well potential. 



By including the charging energies due to the capacitance of the junctions, the Hamiltonian of our circuit 

p2 p2 

H = —f- + —2- + Ej [2 + a - 2 cos(<p p ) cos(^ m ) - a cos(2tt/ + 2<p m )} (5) 

ZlVlp ZIVLui 

where ip p = ipi + ip 2 , ¥>m = <pi - <f2, P P = M p d(p p /dt, P m = M m dip m /dt, M p = (^ /2n) 2 2C, and M m = 
(<i> /27r) 2 2C(l + 2a). The number of degrees of freedom in the problem was reduced by the fluxoid quantization 
condition^ 

2tt$ , . 

<pi + <-Pz - <^2 = 2nn + (6) 

which forces the sum of the gauge invariant phases to be proportional to the amount of flux quanta modulo an integer 
multiple of 2tx. 

We have chosen to associate the capacitive energy, the first two terms in (J5jl , with the kinetic energy, and the ideal 
Josephson energy, the last four terms in JSJ, with the potential energy. The potential energy is that of an infinite 
lattice of double wells, as seen in figure Efb) . The arrow in the plot shows the direction one would take to traverse 
from one side of a double well to another. The barrier between the left and right sides of a single double well is much 
lower than any barrier between different double wells. 

Though quantum mechanics plays a foremost role in deriving the constitutive relations for the superconducting 
circuit elements, the Hamiltonian for the circuit itself so far has been classical. The quantum version of the circuit 
can be understood by imagining a phase "particle" in the potential shown in figure |2{b). The behavior of this 
"particle" is analogous to a particle with an anisotropic mass moving in a 2D periodic potential, and so there exist 
energy bands in a fc-space, which is here related to the charge stored capacitively by the Josephson junctions. By 
properly choosing Ej/E c , where E c — e 2 /2C, one can remove any k (and hence charge) dependence in the energy 
of the system, and hence can reduce the problem to that of an effective double well. What we have done is choose 
parameters such that tunneling between adjacent double wells can be neglected relative to the tunneling within a 
double well in the tight-binding solution, making the solution effectively that of a single double well. 

By considering only the lowest two levels of the double well, the equivalent Hamiltonian is 

H = Q Ip(J- -)a z - ra x (7) 

where ±J p are the eigenvalues of circulating current for the two a z eigenstates and r is the tunneling element from 
one side of the double well to the other. The energies of the two eigenstates along with a sketch of the double well 
as a function of applied flux are shown in figure One significant difference between this qubit and the one used in 
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FIG. 3: The energy levels of the PC Qubit are shown as a function of /. The eigenstates of the system change with / and 
are labeled on the plot. The change in the potential of the phase particle is also depicted at the top of the plot. The energy 
difference between the states at / = 1/2 is seen to be twice the intra- well tunnelling. 

NMRQC is the presence of the a x term. The implication of such a term is that the energy of the eigenstates as well as 
the eigenstates themselves change as the bias field is modified. In figure|j|]we see that at the classical degeneracy point 
/ = 1/2 the qubit's eigenstates are a x eigenstates, while far from / = 1/2, but still far from / = 1, the eigenstates are 
those of a z . The same thing happens for / < 1/2, but now the eigenstates have switched energies, i.e., the ground 
state here is the first excited state of a z and vice versa. 

The PC Qubit has some advantages over other superconducting qubits. Charge fluctuations, a consequence of 
trapped substrate charge, are deemed inconsequential through the choice of parameters used when designing the PC 
Qubit circuit. Also, flux noise has been reduced in this system over other flux qubits since this system has a smaller 
loop. 

A typical conceptual misconception can be addressed at this point. The two different states used in computation are 
not related to single Cooper pair behavior. Rather, they are macroscopically distinct states described by the circulating 
current due to millions of Cooper pairs, characterized by different average induced fluxes when in a magnetic held. 

As seen in section^ the qubits will need to be coupled. For the PC Qubit, just as microwaves can only be coupled 
in through a z , coupling between qubits can only be of the form a z a z . Other coupling terms can be introduced by 
design, even for these planar devices^. 



We now show how one can use the PC Qubit to simulate the ID diffusion equation. In section llV Al we elaborate 
on a scheme based upon changing the flux bias points of the qubits during the algorithm, which will lead to a very 
general initialization scheme, but a less general collision. In section lTVBI we discuss a more general collision, analogous 
to that done in NMRQC, and how to initialize the qubits before this general collision. 



The first of the four steps of the algorithm is initializing each qubit at each node. As discussed in section [HI each 
qubit must be initialized into a state of real and positive phase in its own Hilbert space. This set of states consists 
of all those lying on the real phase geodesic between the ground and first excited states on the Bloch sphere. The 
ground state of the PC Qubit as a function of applied flux coincidentally also occupies exactly this geodesic on the 
Bloch sphere, as discussed in section ITTT1 Initialization can thus be accomplished while staying in the ground state by 
adiabatically changing the applied magnetic flux, as depicted in figure [21 



IV. IMPLEMENTATION WITH THE PC QUBIT 



A. The Multiple Bias Point Implementation 
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Flux Seen By Qubit 2,/ 2 (<I> ) 

FIG. 4: The overlap between the first (second) excited state |1) (|2)) of the PC Qubit coupled system and the 1 01) ( 1 10) ) 
computational state are plotted when qubit 1 is biased at / = 0.508. The coupling between |1) and |2) in the presence of an 
AC magnetic field is also plotted. Qubit coupling equal to r (same for both qubits) was assumed in the calculation. 

The flux used to set the state of one qubit will be affected by the state of the other qubit and its bias current. 
This permanent inductive coupling can be accounted for by slightly adjusting the applied flux to compensate for the 
flux introduced by the other qubit and its bias line. All errors due to approximations made when initializing by 
rotating dynamically on the Bloch sphere, including errors due to decoupling, have been avoided^. We emphasize 
that the initialization portion of the algorithm is identical for any simulation, whether it be for a different equation, 
a multi-phase simulation, or in a different number of dimensions. 

The second step of the algorithm is the collision. Here we study a very specific unitary transformation, the y/swap 
described in section m This matrix simply "half-way" swaps the middle two (first and second excited) computational 
states of the coupled system. In NMRQC, the coupled eigenstates are exactly those computational states, but there 
are no direct matrix elements connecting these states^. When the PC Qubits are coupled, the first and second excited 
states of the four- level system, denoted as |1) and |2) respectively, are in general not the same as the computational 
basis states the y/swap intends to affect. However, the dc bias fields of each qubit can be tuned to make these two 
sets of eigenstates coincide. Once this is done, one can then implement the ^Jswap by simply oscillating the magnetic 
field bias at the frequency corresponding to the energy difference between the middle two eigenstates. This is just 
a Rabi oscillation between the middle two eigenstates, and since one wants to only "half-way" swap the states, the 
radiation should only be left on for a quarter of a Rabi period. 

Besides finding the appropriate bias points such that the middle two eigenstates of the coupled system are very 
similar to the middle two computational states, one must also verify that the coupling between these states in the 
presence of an oscillating magnetic field is non-zero. The results of these calculations are shown in figure 0] The bias 
point of qubit 1, fi, must be chosen to be far from 1/2, but not too far. In these calculations we take /i = 0.508. In 
the figure we see that when qubit 2 is biased at around f% = 0.51, the first two system excited states are very similar 
to the middle two computational states, with overlap elements of about 0.97. At this same bias point one sees a Rabi 
matrix element of about 0.02, which is more than sufficient for our purposes. 

This approximate swap has been incorporated into simulation of the FQLGA for the ID diffusion equation and 
the results are pictured in figure [SJ Snapshots of three different times have been shown, for both an ideal simulation 
and one including the error introduced due to the approximate collision. At time t = one can see that we have 
initialized our fluid to a gaussian profile. Later time steps of the ideal implementation show the expected spreading 
due to diffusion, while conserving the total number of particles. Increase in the diffusion constant of the approximate 
collision when compared to the ideal simulation results from the enhanced population in the |00), 1 01) , and 1 10) states 
relative to the 1 1 1) state due to extra matrix elements in the approximate swap that couple the four states. 

Even with an ideal swap operator, an interesting timing issue arises upon non-adiabatically switching the bias 
fields from the initialization settings to the proper settings to do a Rabi oscillation between the two middle product 
states. We first illustrate this timing issue and then show how it can be made negligible by making a larger ensemble 
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FIG. 5: The results of the FQLGA are simulated having accounted for the approximate nature of the collision proposed in 
section II V Af +). Qubit coupling equal to r (same for both qubits) was assumed in the approximate swap simulation. The ideal 
results of the FQLGA are also shown(o). 



measurement. 

Once the applied fluxes are changed to those appropriate to perform the approximate swap, the initialized states 
will most likely not be eigenstates anymore, and hence will begin to precess due to a time-independent perturbation. 
Assuming things can not be accurately controlled at these timcscalcs, one will have now introduced a random phase 
difference between the two qubits due to this Larmor precession. This effect is pictured in figure EJ The states before 
the bias fields are switched lie along the same geodesic. Upon changing the magnetic flux seen by each qubit, the 
qubits begin to precess, out of phase. 

The effect of this phase difference 6 on the algorithm will be to alter the fraction of particles at each lattice site, 
post-collision, that are "moving" to the right and to the left. The results of measuring the post-collision occupation 
probabilities having accounted for a constant phase difference is summarized by 

Pi =Pi. s=0+jsm(5) (8a) 
P2 = P 2 . 4=o-7sin((S). (8b) 

The effect of this error on the simulation is effectively averaged away when an ensemble is measured, since 5 is randomly 
different for each member of the ensemble. These results are shown in figure [7] One can see small random deviations 
from the ideal simulation that can be made infinitesimally small by measuring a larger ensemble (an ensemble average 
of 1000 repeated measurements was simulated here). 

In summary, an initialization scheme has been developed that is not available to qubits with only one term in 
their Hamiltonian. This initialization scheme is limited only by the precision of the current source used to create 
the magnetic field that biases the qubit. The scalability of this scheme relative to those used in NMRQC is an 
interesting question, but is not resolved here. The collision implementation is also unique to qubits with multiple 
term Hamiltonians, but the unitary transform implemented is unique to the diffusion equation, and fortuitously 
simple. A collision scheme that could be generalized to any unitary transformation would be much more useful. 



B. Generalized NMRQC-like Implementation 

Generalization of the above implementation to any fluid dynamics, i.e., any unitary transformation, can be done in 
an analogous way to NMRQC schemes. Generalization of the collision transformation consists of using a universal set 
of quantum computation gates, and decomposing all transformations into a sequence of these 2 . In NMRQC collision 
is performed by a sequence of single qubit unitary transformations and coupled free evolution. In this section we will 
begin by discussing the single qubit rotations needed for a general decomposition, and briefly mention the role they 
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FIG. 6: The unfilled circle and triangle represent two typical initialized states at one lattice point, both on the same north 
pole to south pole geodesic, before their flux bias is changed to perform the collision. The filled circle and triangle represent 
the same states after imprecise bias changing has occurred. Imprecisely timed Larmor precession introduces a random phase 
difference S between the two states. 
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FIG. 7: The results of the FQLGA are simulated having accounted for a random phase difference introduced before the collision 
for each member of the measurement ensemble(+). The ideal results of the FQLGA are also shown(o). 



could play in initialization. We will then explore the free evolution of a coupled PC Qubit system, and then show 
how to combine the single and coupled pulses to implement the collision of the ID FQLGA for the diffusion equation. 

Single qubit transformations can most easily be achieved in a rotating frame, since here the frequency of precession 
can be much lower than the Larmor time scale. For this implementation we will only study the case where our qubit 
is biased at / = 1/2. This discussion is easily generalized to any bias point, but the mathematical notation can get 
quite cumbersome. The Hamiltonian of the PC Qubit in an applied AC field is 

H = w I x + g a cos(w t + 4>)I Z (9) 

where w is the frequency of the applied field, g Q is proportional to the amplitude of the applied field, 4> is the phase 
of the applied field, and li — hbijl. 
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FIG. 8: The results of the FQLGA are simulated for NMR-like single qubit pulse initialization, where errors arise from 
initializing a coupled ground state that is not a product state(+). The ideal results of the FQLGA are also shown(o). 



In the frame rotating about I x this Hamiltonian becomes 

H = ~g [cos{4>)I z + sin(0)4] . (10) 

The quantum state will now precess in this frame about the axis defined by 0, with the angle through which the state 
has precessed given by 9 = g t/2. It will be convenient to only consider the set of two rotations defined by <f> = and 
7r/2, which are rotations about I z and I y respectively. These rotations, denoted as R z (0) and R y {9) respectively, can 
be used in conjunction to bring the qubit state to any point on the Bloch sphere in the rotating frame. 

One can use these single qubit rotations not only as part of the collision, but also for initializing, since they can 
bring the qubit state to anywhere on the Bloch sphere. As already discussed in section IIV Al the ground state of 
a coupled PC Qubit system is not in general the product of single qubit ground states. Thus, when initializing a 
qubit via R z {6) and R y {6), one is not starting rotation from the single qubit ground state. However, since the ground 
state is very close to a product of single qubit ground states, this difference is nearly negligible. In figure we show 
the effects of incorporating this error into the algorithm when the coupling constant is taken to be 1/10 of the qubit 
resonant frequencies, a rather exaggerated estimate since the coupling is usually much smaller. The diffusion constant 
is decreased by this approximation, due to the enhanced population in the |11) state relative to the |00) state from 
the coupling. 

The other gate operation needed to form a universal set for general decomposition is coupled free evolution. The 
first thing to do is go to a co-rotating frame where one can have a coupled Hamiltonian only, i.e., no single qubit 
terms, that is time- independent. In NMRQC this is done by going to the frame where both qubits are rotated around 
the z axis. However, since our coupling does not commute with our single qubit terms, a different method will be 
used. For notational convenience only, we consider the case where both qubits are biased at / = 1/2, where our 
Hamiltonian is 

h = wlil + w 2 Jl + y Ji 2 [iii 2 2 ]. (ii) 

In the co-rotating frame where both qubits are rotating around the x axis, one has the Hamiltonian 

H = ^j 12 [ili 2 z +i^} (12) 

as long as wl = w 2 . This constraint of w\ — w 2 imposes limitations on some NMRQC initialization schemes which 
use frequency selective initialization. 

One can now rewrite the unitary collision transformation in the following suggestive way: 

yfswap = exp[-i^(<7^ + c^)] exp[-^(T;^] (13) 
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The first term is just free evolution in the co-rotating frame. The second term can be written as: 

expH|<7^] = Riq)Rl(l)exp[-i^alal]Rl(~)Rl(~) (14) 

where the middle term can be written as: 

cxpH^a 2 2 ] =expH^a 2 2 + ^)]i?i(^)expH^ (15) 

Hence one can perform a decomposition of the collision transformation into a sequence of single qubit rotations and 
coupled free evolution. 

In summary, we have shown that the PC Qubit can implement the unitary transform that performs collisions in 
the ID FQLGA for the diffusion equation by a single and coupled qubit evolution decomposition. The single qubit 
rotations were shown to be feasible for qubit initialization as well, with a slight approximation due to the coupled 
ground state that is not a product state. The coupled free evolution was seen to require identical qubit frequencies 
over a lattice site, making initialization a bit more challenging. 



V. CONCLUSIONS 



In this paper we have shown that the implementation of the FQLGA for the ID diffusion equation is feasible 
with PC Qubits. We began by considering the simplest scheme possible using the PC Qubit. This consisted of first 
initializing the qubits while keeping them in their ground state, and then performing the collision by quickly changing 
their flux bias points and then performing a single n/2 pulse. This initialization technique could prove useful, but 
the way we have implemented the collision is not easily generalized to other collisions. We needed to develop a more 
general collision scheme, and then see how we could initialize in conjunction with that new scheme. 

A more general collision transformation was then discussed by decomposing the unitary matrix into a sequence of 
single qubit rotations and coupled free evolution. We first developed single qubit rotations for the PC Qubit that could 
be used as part of the collision decomposition as well as for initializing the occupation probabilities. The initialization 
was considered only approximate due to the permanent non-commuting coupling between qubits. For the coupled free 
evolution we saw that transforming to a rotating frame analogously to NMRQC set a strong but feasible constraint 
on the frequencies of our qubits. Ultimately one would like to remove the constraint of equal frequencies, so that 
frequency-selective initialization can be done analogously to the NMRQC implementation, alongside the very general 
collision scheme. One would then also need to account for initialization pulses rotating states from a non-product 
ground state. 
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